Uniaxial Strain Dependence on Angle-Resolved Optical Second Harmonic Generation from a Few Layers of Indium Selenide

Indium selenide (InSe) is an emerging van der Waals material, which exhibits the potential to serve in excellent electronic and optoelectronic devices. One of the advantages of layered materials is their application to flexible devices. How strain alters the electronic and optical properties is, thus, an important issue. In this work, we experimentally measured the strain dependence on the angle-resolved second harmonic generation (SHG) pattern of a few layers of InSe. We used the exfoliation method to fabricate InSe flakes and measured the SHG images of the flakes with different azimuthal angles. We found the SHG intensity of InSe decreased, while the compressive strain increased. Through first–principles electronic structure calculations, we investigated the strain dependence on SHG susceptibilities and the corresponding angle-resolved SHG pattern. The experimental data could be fitted well by the calculated results using only a fitting parameter. The demonstrated method based on first–principles in this work can be used to quantitatively model the strain-induced angle-resolved SHG patterns in 2D materials. Our obtained results are very useful for the exploration of the physical properties of flexible devices based on 2D materials.


Introduction
Two-dimensional (2D) materials have drawn much attention since the discovery of graphene in 2004 [1]. Subsequently, graphene and the subsequent derivative 2D materials have been demonstrated for many applications such as in electronic and optoelectronic devices [2,3]. Since van der Waals materials can be mechanically exfoliated and transferred to arbitrary substrates, including bendable polymeric or plastic substrates, studies of the strain-dependent properties of 2D materials are therefore beneficial for the development of flexible devices [4,5]. Indium selenide (InSe) is an emerging star in van der Waals semiconductors [6] due to its superior properties such as ultrahigh mobility [7,8] and large elastic deformability [9,10]. It has been applied to bendable photodetectors with a good performance and a broad spectral response [11]. For the application of flexible devices, it is important to understand how strain affects the relevant physical properties. For example, InSe was found to exhibit the ultrasensitive tunability of the direct bandgap by applying uniaxial strain [12,13]. The variation in optical properties under strain could be utilized as a tool to map the intensity distribution of strain such as photoluminescence spectroscopy and micro-Raman spectroscopy [14,15]. Recently, optical second harmonic generation (SHG) techniques have been demonstrated to be a powerful and unique tool for mapping strain vectors. In addition to the intensity distribution of strain, the direction of the strain distribution can also be mapped by analyzing the angle-resolved SHG patterns [16].
SHG is a second-order nonlinear process that only occurs in crystals without inversion symmetry. In a thick crystal, the conversion efficiency of SHG is maximized in the condition of phase matching. For thin crystals with a thickness below 100 nm, a few methods have been proposed to increase the efficiency of SHG recently [17][18][19]. For example, an optical interference layer at the bottom of the thin crystals enhances the SHG when the thickness of the interference layer is designed in the phase condition for the constructive interference of both the fundamental and second harmonic frequencies [17]. Nanostructured crystals serving as dielectric nanoantennas for second harmonic emissions have also been demonstrated [18,19]. While the Mie resonant modes are matched with the fundamental or second harmonic frequencies, the efficiency of SHG is increased. The nonlinear optical process highly depends on the motion of the electrons in the lattice structure. It is thus capable of determining the orientations of single crystals [20]. For example, monolayer transition metal dichalcogenide (TMD) materials, such as MoS 2 , WS 2 , and WSe 2 possess hexagonal structures with a point group of D 3h . When a linearly polarized laser is normally incident to a TMD monolayer, SHG light with the same polarization can be recorded as a function of the rotation angles of the TMD. This angle-resolved SHG pattern exhibits six-fold rotational symmetry, and the maximum SHG intensity occurs during laser polarization along the arm-chair (AC) direction [21][22][23]. If strain is applied to the hexagonal TMD monolayers, the crystal symmetry reduces, and the angle-resolved SHG patterns vary accordingly [16,24,25].
Since not only the intensity of SHG, but also the symmetry of the angle-resolved SHG patterns are analyzed, it is possible to monitor the local strain vector. InSe possesses high mobility, which is an advantage over the TMDs in certain device applications. It possesses a periodic honeycomb lattice in three specific phases (β, γ, and ε) based on different stacking sequences [26]. SHG is allowed in monolayers and a few layers of InSe in γ phase/ε phase due to their non-centrosymmetric structures [27][28][29][30].
In this work, we investigated the strain dependence on angle-resolved SHG from mechanically exfoliated γ-InSe flakes. We mapped the SHG images of the strained InSe flakes with different rotation angles and studied the angle-resolved SHG pattern. Overall, the SHG intensity decreased with increasing compressive strain. By performing firstprinciples electronic structure calculations within the density functional theory, we also studied the uniaxial strain dependence on SHG susceptibilities. The calculated susceptibility decreases with increasing compressive strain, which agrees with the experimental results. We demonstrate a method based on first-principles calculations to quantitatively model the strain-induced angle-resolved SHG patterns in InSe flakes, which is applicable to other 2D materials. This model is useful to analyze angle-resolved SHG patterns for the strain mapping of 2D materials.

Growth of InSe Single Crystals and Their Characterizations
InSe crystals were grown by the vertical Bridgman method. We purchased 99.999% pure Indium (In) and Selenide (Se) chemicals from Sigma Aldrich. The growth of the InSe single crystals was performed in quartz ampoules by mixing In and Se powders under a vacuum of 10 −4 Pa. The ampoules were placed in a horizontal furnace at 550 • C for two days to obtain a partially mixed crystal. Later, the mixed crystals were used to grow the pure InSe single crystal with the help of a vertical Bridgman setup. The furnace during the Bridgman process was maintained at a high temperature of about 850 • C. The mixed crystal in the ampoule was hung and heat treated at 850 • C for 24 h; once the mixed crystal melted entirely in the ampoule, the ampoules were lowered through a temperature gradient of 1 • C at a rate of 0.1 mm/h to obtain a large-sized single crystal. Using this procedure, we obtained over 5 cm long and 1 cm in diameter InSe single crystals in bulk.
Powder X-ray diffraction (PXRD) was employed to study the lattice phases. The result of the Rietveld refinement of X-ray diffraction was carried out by the software GSAS, while the high orientation of the powdered InSe was taken into account for the refinement. The lattice structure was determined as γ-InSe (R3m) with lattice constants a = b = 4.006 Å, and c = 24.969 Å, respectively. The crystal phase was also confirmed by oblique incident SHG measurements [31]. Although SHG is allowed in both the γ phase and ε phase due to their non-centrosymmetric structures, their different stacking orders lead to distinguishable zrelated SHG susceptibilities, which are measurable through the angle-resolved p-polarized SHG (p-out) with an oblique incidence of the p-polarized fundamental light (p-in). The p-out p-in SHG angle-resolved pattern of our InSe crystal manifested 3-fold rotational symmetry, indicating the feature of γ phase because the p-out p-in angle-resolved SHG pattern of ε phase InSe exhibits 6-fold rotational symmetry.
The sample structure of γ-InSe is shown in Figure S1 in the Supplementary Material. The cross-section view along the [010] direction in Figure S1a indicates how the monolayers of InSe are stacked in an ABC manner. The top view of γ-InSe along the [001] direction is shown in Figure S1b, where the AC and zigzag (ZZ) directions are labelled with black arrows. The bulk γ-InSe crystal belongs to the C 3v point group and exhibits a non-centrosymmetric structure.

Sample Preparation and Strain Engineering
Flakes of InSe in a few layers were prepared by the mechanical exfoliation of the γ-InSe crystal with polydimethylsiloxane (PDMS). Multiple flakes were subsequently transferred to a polyethylene naphthalate (PEN) flexible substrate, with an SU-8 photoresist film on top. The thickness of the PEN substrate was about 125 µm, and the thickness of SU-8 photoresist film was in the order of 1 µm. On top of the PEN substrate with the SU-8 photoresist film, the InSe flakes were examined using an atomic force microscope (AFM). Figure S2a in the Supplementary Material shows the AFM image around the edge of an InSe flake. Four line profiles of height were analyzed, as shown in Figure S2b, which reveal that the thickness is in the range between 7.3 and 8.7 nm. The thickness of the InSe flake we studied in this work is~8 nm, which corresponds to~10 layers.
We referred to the experimental configuration in the literature [32] to design a device for applying uniaxial strain. As shown in Figure S3 in the Supplementary Material, the strain device was combined with two sliding stages and screws. The flexible PEN substrate was placed between the two sliding stages. Uniaxial strain was applied by curving the PEN substrate via spinning the screws to move the two sliding stages toward each other. Figure S4 in the Supplementary Material shows the schematic of the curved substrate under uniaxial strain. The uniaxial strain ε in this work was determined following the relationship ε = d/2R, where d is the thickness of the substrate and R is the radius of curvature of the substrate. In the experiments, 125 µm was used for d. To obtain the radius of curvature, the cross-section images of the curved PEN substrate on the strain device were recorded using a digital camera. The error of fitting the radius of curvature and retrieving the strain was within 3%.

SHG Measurements
The schematic of the experimental setup is shown in Figure S5 in the Supplementary Material. The ultrafast pulses centered at~790 nm were generated by a Ti:sapphire laser at 80 MHz repetition rate (Newport, Tsunami). The bandwidth of the pulses was~12 nm. A set of a half-wave plate (HWP) and a polarizing beam splitter (PBS) cube was used to control the power to prevent damage occurring to the samples. A pulse compressor was used to compensate the dispersion of the transparent optical components to maximize the efficiency of SHG from the samples. The beam was directed through a galvo-mirror system and a long working distance objective lens (20×, NA = 0.35) to scan the samples in two dimensions. A linear polarizer (LP) was used to purify the polarization of the incident beam. A dichroic mirror was used to separate the fundamental incident beam and SHG beam from the sample. The SHG photons were sensed by a single-photon-counting photomultiplier tube (PMT). A 400 nm bandpass filter (BPF) was placed in front of the PMT to eliminate the light at the fundamental frequency. A linear polarizer was placed in front of the bandpass filter, which was rotated to selectively measure the SHG light with polarization that was parallel or perpendicular to the polarization of the fundamental incident light. A motorized rotational stage (Newport, URS150BCC) was on an inverted microscope for recording the angle-resolved SHG images with a stepping angle of 10 degrees.

Theoretical Calculations
To theoretically obtain the geometrical structures of InSe under strain and calculate the SHG susceptibilities of the strained InSe, we performed first-principles calculations based on the density functional theory (DFT) using the OpenMX code [33], where generalized gradient approximations (GGA) [34], norm-conserving pseudopotentials [35], and optimized pseudoatomic basis functions [36] were adopted. We constructed a lattice of the γ phase InSe, with a unit cell consisting of three InSe monolayers in an ABC stacking order, as shown in Figure S1 in the Supplementary Material. The x-axis was set to be the AC direction. For the atomic basis functions, three, two, and two optimized radial functions were allocated to the s, p, and d orbitals, respectively. We adopted a cutoff radius of 7 bohr for the In atom, denoted as In7.0-s3p2d2. Similarly, for the Se atom, Se7.0-s3p2d2 was adopted. A cutoff energy of 500 Ry was used for numerical integrations and for the solution of the Poisson equation. An 8 × 8 × 2 k-point sampling was used. The structures were fully relaxed without considering spin-orbit coupling at the specified strain, that is, only the components of the lattice vectors associated with the desired strain were fixed during the calculations. For the unstrained ground-state structure, no constraint was applied to the lattice vectors. Once the optimized lattice constants were obtained, the atomic positions were relaxed again with the effect of spin-orbit coupling. All of the atomic forces in the resulting lattice satisfied the criterion of being smaller than 0.0001 Ha/bohr. The electronic self-consistent field (SCF) criterion was set to 10 −9 Ha, which guaranteed that the difference between the band energies in the last two consecutive SCF steps was smaller than the criterion.
The spectra of SHG susceptibilities were calculated via the momentum matrix elements using the basis of the adopted pseudoatomic basis functions [37], following the formula derived by Sipe et al. for the independent-particle approximation [38,39]. A 60 × 60 × 8 k-point sampling was used for obtaining the SHG susceptibilities with a broadening parameter of η = 0.1 eV. The electronic temperature was set to 300 K for all of the SCF and SHG calculations based on the Fermi-Dirac distribution function. Due to the band gap error in the DFT calculations [40], the calculated band gaps were underestimated, and a scissor operator with 0.6487 eV was applied to increase the band gaps by shifting the eigenvalues of unoccupied Kohn-Sham orbitals. The modification to the momentum matrix elements, which was accompanied by the energy shift, was also taken into account in the SHG calculations [39].  Figure 1 were SHG signals from the InSe flakes, we measured the spectrum of the collected light in front of the PMT, as shown in Figure S6a in the Supplementary Material. The Gaussian-like feature centered at 394 nm confirmed the SHG of the incident light pulse centered at~790 nm. We also conducted power-dependent measurements to confirm the SHG process, as shown in Figure S6b in the Supplementary Material. The fitting line with a slope of~2 in the log-log scale confirmed the quadratic dependence of SHG intensity on the fundamental laser power. For the SHG images in Figure 1, the polarizers were set so that the polarization of the detected SHG was parallel to the polarization of the incident light. While the samples were rotated from −20 to 30 degrees, the optical intensity of the biggest flakes among the three flakes increased to a maximum at 0 degree, shown in Figure 1c, and decreased to a minimum at 30 degrees, shown in Figure 1f. To further quantitatively analyze the angle-resolved SHG intensity, we summed up the signals over the biggest flakes in the images of different rotational angles. The angular dependence on the SHG intensity of parallel polarization (denoted as I (2ω) ) from the biggest InSe flake is revealed by the blue squares in Figure 2. Because the hexagonal γ-InSe exhibits the point group of C 3v , the angle-resolved I (2ω) under normal incidence also reveals six-fold rotational symmetry. The maximum intensity occurs at the angle of the incident light polarization along the AC direction [41]. The angle-resolved I (2ω) was fitted with cos 2 3θ, as shown by the blue line in Figure 2. Here, θ is defined as the angle of AC direction with respect to the linear polarization of the incident laser beam. Angle-resolved analysis is capable of determining the crystal orientations of different flakes on the substrate. For example, while the SHG intensity of the upper left flake in Figure 1a was traced for different rotational angles in Figure 1a-f, the maximum SHG intensity occurs at 10 degrees in Figure 1d. In contrast to the maximum SHG intensity of the biggest flake occurring at 0 degree in Figure 1c, the orientation of the upper left flake was rotated by 10 degrees with respect to that of the biggest flake. We also conducted mapping measurements with SHG polarization perpendicular to the polarization of the incident light (denoted as I

Results and Discussion
. For a similar angle-resolved analysis of the biggest flake, the orange points in Figure 2 depict the experimental data in the polar plot of I (2ω) ⊥ , which manifest a thirty-degree-rotated pattern of I (2ω) , and they are fitted with sin 2 3θ.
ization of the incident light. While the samples were rotated from −20 to 30 degrees, the optical intensity of the biggest flakes among the three flakes increased to a maximum at 0 degree, shown in Figure 1c, and decreased to a minimum at 30 degrees, shown in Figure  1f. To further quantitatively analyze the angle-resolved SHG intensity, we summed up the signals over the biggest flakes in the images of different rotational angles. The angular dependence on the SHG intensity of parallel polarization (denoted as ∥ ( ) ) from the biggest InSe flake is revealed by the blue squares in Figure 2. Because the hexagonal -InSe exhibits the point group of , the angle-resolved ∥ ( ) under normal incidence also reveals six-fold rotational symmetry. The maximum intensity occurs at the angle of the incident light polarization along the AC direction [41]. The angle-resolved ∥ ( ) was fitted with cos 3 , as shown by the blue line in Figure 2. Here, is defined as the angle of AC direction with respect to the linear polarization of the incident laser beam. Angle-resolved analysis is capable of determining the crystal orientations of different flakes on the substrate. For example, while the SHG intensity of the upper left flake in Figure 1a was traced for different rotational angles in Figures 1a-f, the maximum SHG intensity occurs at 10 degrees in Figure 1d. In contrast to the maximum SHG intensity of the biggest flake occurring at 0 degree in Figure 1c, the orientation of the upper left flake was rotated by 10 degrees with respect to that of the biggest flake. We also conducted mapping measurements with SHG polarization perpendicular to the polarization of the incident light (denoted as ( ) ). For a similar angle-resolved analysis of the biggest flake, the orange points in Figure 2 depict the experimental data in the polar plot of ( ) , which manifest a thirtydegree-rotated pattern of ∥ ( ) , and they are fitted with sin 3 .    Figure 3. The best fitting lines with cos 2 3θ are shown to guide the eyes. Overall, the intensity of SHG decreases with increasing compressive strain.  In order to obtain more insight into strain-induced angle-resolved SHG patterns, we theoretically calculated the uniaxial strain dependence on SHG susceptibilities. First, we used the first-principles method to calculate the band structures of InSe under different compressive strains. The strain was applied along the AC direction, which is defined as the x direction. Subsequently, the spectra of SHG susceptibilities were calculated following the method described in Section 2.4 and are shown in Figure S7 in the Supplementary Materials. The SHG susceptibilities are dependent on the photon energy. Here, we focus on the non-zero SHG susceptibilities at 1.55 eV, which is the photon energy of the incident light in the experiments. Figure 4a shows the strain dependence on ( ) at 1.55 eV. The real parts, imaginary parts, and phases of non-zero SHG susceptibilities at 1.55 eV are shown in Figure S8 in the Supplementary Material. When the compressive strain increases from 0% to −2%, all of the non-zero components of ( ) decrease. These trends agree with the experimental results that the SHG intensity decreases with increasing strain. We developed a model based on these first-principles results to quantify the rotational angle In order to obtain more insight into strain-induced angle-resolved SHG patterns, we theoretically calculated the uniaxial strain dependence on SHG susceptibilities. First, we used the first-principles method to calculate the band structures of InSe under different compressive strains. The strain was applied along the AC direction, which is defined as the x direction. Subsequently, the spectra of SHG susceptibilities were calculated following the method described in Section 2.4 and are shown in Figure S7 in the Supplementary Materials. The SHG susceptibilities are dependent on the photon energy. Here, we focus on the non-zero SHG susceptibilities at 1.55 eV, which is the photon energy of the incident light in the experiments. Figure 4a shows the strain dependence on χ (2) at 1.55 eV. The real parts, imaginary parts, and phases of non-zero SHG susceptibilities at 1.55 eV are shown in Figure S8 in the Supplementary Material. When the compressive strain increases from 0% to −2%, all of the non-zero components of χ (2) decrease. These trends agree with the experimental results that the SHG intensity decreases with increasing strain. We developed a model based on these first-principles results to quantify the rotational angle dependence on the polarized SHG intensity and to compare them with the experimental data. The second-order polarization P (2ω) can be expressed as (1) ijk is the second-order susceptibility tensor component and E j (E k ) is the electric field. The indices i, j, and k can be 1, 2, and 3 corresponding to x, y, and z in Cartesian coordinates, respectively. Because E j and E k have the same frequencies ω in SHG, the second-order susceptibility possesses intrinsic permutation symmetry, which is χ (2) ijk = χ (2) ikj . By considering this permutation rule, the third-rank tensor χ (2) for SHG is contracted to a 3 × 6 matrix with element d im . The index m can be 1 = xx, 2 = yy, 3 = zz, 4 = yz, zy, 5 = xz, zx, and 6 = xy, Nanomaterials 2023, 13, 750 7 of 11 yx. Since the unstrained γ phase InSe belongs to the C 3v point group, its d-matrix has eight non-zero components, as follows [20].
where −d 11 = d 12 = d 26 , d 15 = d 24 , and d 31 = d 32 . Because the fundamental beam is normally incident to InSe, we consider the polarization of light on the x-y plane with θ respect to the x-axis. The electric field vector can be written as [20] (3) Nanomaterials 2023, 13, x FOR PEER REVIEW 8 of 11 without strain, the overall intensity reduces, but it exhibits anisotropy due to the strain. The rotational symmetry is reduced from six-fold to two-fold symmetry because the uniaxial compressive strain along the AC direction lowers the crystal symmetry of γ-InSe from the point group of to the point group of . The two petals at 0° and 180°, which mean that the laser polarization is along the uniaxial strain direction, reduce more than the other four petals do. In Equation (6), the cos 3 term, which is a function with six-fold rotational symmetry, describes the angle-resolved ∥ ( ) without strain, as shown by the blue line in Figure 4b. The modification due to the strain effect is found in the term cos + ( + 2 ) sin cos , which is a function with two-fold rotational symmetry.
Therefore, the modification due to the strain term reduces the symmetry of the angleresolved ∥ ( ) , as shown by the red line in Figure 4b. To quantitatively understand the appropriateness of the aforementioned model, we used the strain values in the experiments to calculate the angle-resolved ∥ ( ) in Equation (6) for a comparison with the experimental data. Note that , , , and are derived from the results of first-principles calculations in Figure 4a and are not the fitting parameters. With only a scaling factor , which is multiplied on the right-hand side of Equation (6), the blue line in Figure 5a reveals the best fit curve for = 0. By using the same scaling value of , the red line depicts ∥ ( ) for = −0.83%. Although the intensity of both the experimental and theoretical data drops due to the strain, that of the experimental data drops more, as shown in Figure 5a. To quantitatively evaluate the deviation of theoretical curves from the experimental results, the r-squared values of each theoretical curve for =−0.18%, −0.36%, −0.53%, and −0.83% were calculated as 0.992, 0.987, 0.985, and 0.968, respectively, as shown in Figure S9a-d in the Supplementary Material. The r-squared value decreases with an increasing strain value, indicating that the theoretical curves deviate more from experimental results under higher strain. In other words, the model based on the first-principles calculation underestimates the strain-induced variation in the SHG susceptibilities.
The model should be modified to appropriately fit the experimental data. Here, we additionally introduced a correction parameter and replaced in Equation (6) with . By fitting the parameter = 4.48, we obtained the maximum summation of the r- The parallel SHG intensity for unstrained InSe is The coefficients A, B, and C are complex numbers and specify the linear deviations of complex SHG susceptibilities d 11 , d 12 , and d 26 with respect to strain, respectively. Since (6) d is the calculated SHG susceptibility without strain. The coefficients A, B, and C were obtained by generating linear fits of the calculated strain dependence on the SHG susceptibilities shown in Figure S8a,b in the Supplementary Material. The blue line and the red line in Figure 4b reveal the calculated angle-resolved I (2ω) for ε = 0 and ε = −20% in Equation (6), respectively. Here, a large strain value of −20% was intentionally used to highlight the variation in the angle-resolved SHG pattern. Compared with the pattern without strain, the overall intensity reduces, but it exhibits anisotropy due to the strain. The rotational symmetry is reduced from six-fold to two-fold symmetry because the uniaxial compressive strain along the AC direction lowers the crystal symmetry of γ-InSe from the point group of C 3v to the point group of C 2v . The two petals at 0 • and 180 • , which mean that the laser polarization is along the uniaxial strain direction, reduce more than the other four petals do. In Equation (6), the d cos 3θ term, which is a function with six-fold rotational symmetry, describes the angle-resolved I (2ω) without strain, as shown by the blue line in Figure 4b. The modification due to the strain effect is found in the term A cos 3 θ + (B + 2C) sin 2 θ cos θ , which is a function with two-fold rotational symmetry. Therefore, the modification due to the strain term reduces the symmetry of the angle-resolved I (2ω) , as shown by the red line in Figure 4b.
To quantitatively understand the appropriateness of the aforementioned model, we used the strain values in the experiments to calculate the angle-resolved I (2ω) in Equation (6) for a comparison with the experimental data. Note that d, A, B, and C are derived from the results of first-principles calculations in Figure 4a and are not the fitting parameters.
With only a scaling factor f 2 , which is multiplied on the right-hand side of Equation (6), the blue line in Figure 5a reveals the best fit curve for ε = 0. By using the same scaling value of f 2 , the red line depicts I (2ω) for ε = −0.83%. Although the intensity of both the experimental and theoretical data drops due to the strain, that of the experimental data drops more, as shown in Figure 5a. To quantitatively evaluate the deviation of theoretical curves from the experimental results, the r-squared values of each theoretical curve for ε = −0.18%, −0.36%, −0.53%, and −0.83% were calculated as 0.992, 0.987, 0.985, and 0.968, respectively, as shown in Figure S9a-d in the Supplementary Material. The r-squared value decreases with an increasing strain value, indicating that the theoretical curves deviate more from experimental results under higher strain. In other words, the model based on the first-principles calculation underestimates the strain-induced variation in the SHG susceptibilities. The model should be modified to appropriately fit the experimental data. Here, we additionally introduced a correction parameter f 2 and replaced ε in Equation (6) with f 2 ε. By fitting the parameter f 2 = 4.48, we obtained the maximum summation of the r-squared values for all of the data with different strains. The r-squared values of each modified theoretical curve for ε = −0.18%, −0.36%, −0.53%, and −0.83% are 0.993, 0.996, 0.996, and 0.995, respectively, as shown in Figure S9e-h in the Supplementary Material. The fitting results for all of the experimental data are also depicted in Figure 5b. By only using a correction parameter, the results from the model based on first-principles calculations quantitatively agree well with the experimental data.

Conclusions
We mapped the SHG images of γ-InSe flakes under uniaxial strain. By recording the SHG intensity mapping of InSe flakes with different azimuthal angles, the angle-resolved SHG patterns were obtained. We also performed first-principles calculations to study the strain dependence on the SHG susceptibilities and the corresponding angle-resolved SHG patterns. Both the experimental and calculated results indicate that the SHG intensity decreased with increasing uniaxial strain on the InSe flake. Our experimental investigation and theoretical model provide a way to analyze angle-resolved SHG patterns for the strain mapping of γ-InSe flakes, which is very useful for better understanding the physical properties of 2D materials implemented in flexible devices. The modeling method proposed in this work is expected to be applicable to investigations of the strain-dependent SHG in other 2D materials.  Figure S3. Schematic of the strain device. Figure S4. Schematic diagram of curved substrate under uniaxial strain. Figure S5. Schematic of the experimental setup. Figure S6. (a) The SHG spectrum with peak wavelength at 394 nm. (b) The power dependent SHG intensity. The slope of the linear fitting curve is around 2 in the log-log scale, indicating a quadratic relationship between the incident power and the SHG intensity. Figure S7. Calculated SHG susceptibility components of γ-InSe under strain of (a) −2% and (b) 0%. Figure S8. Strain dependence on the (a) real part, (b) imaginary part, and (c) phase of calculated , , and at 1.55 eV. Figure S9. R-squared values of the strain-dependent experimental data by using the model from the first-principles methods (a)-(d) without parameters and (e)-(h) with one parameter to linearly modify the strain.

Conclusions
We mapped the SHG images of γ-InSe flakes under uniaxial strain. By recording the SHG intensity mapping of InSe flakes with different azimuthal angles, the angle-resolved SHG patterns were obtained. We also performed first-principles calculations to study the strain dependence on the SHG susceptibilities and the corresponding angle-resolved SHG patterns. Both the experimental and calculated results indicate that the SHG intensity decreased with increasing uniaxial strain on the InSe flake. Our experimental investigation and theoretical model provide a way to analyze angle-resolved SHG patterns for the strain mapping of γ-InSe flakes, which is very useful for better understanding the physical properties of 2D materials implemented in flexible devices. The modeling method proposed in this work is expected to be applicable to investigations of the strain-dependent SHG in other 2D materials.  Figure S3. Schematic of the strain device. Figure S4. Schematic diagram of curved substrate under uniaxial strain. Figure S5. Schematic of the experimental setup. Figure S6. (a) The SHG spectrum with peak wavelength at 394 nm. (b) The power dependent SHG intensity. The slope of the linear fitting curve is around 2 in the log-log scale, indicating a quadratic relationship between the incident power and the SHG intensity. Figure S7. Calculated SHG susceptibility components of γ-InSe under strain of (a) −2% and (b) 0%. Figure S8. Strain dependence on the (a) real part, (b) imaginary part, and (c) phase of calculated d 11 , d 12 , and d 26 at 1.55 eV. Figure S9. R-squared values of the strain-dependent experimental data by using the model from the first-principles methods (a-d) without parameters and (e-h) with one parameter to linearly modify the strain.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.